import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import random

# 生成a到b之间的随机数
def GetRandomFactor():
    random.seed()  # 使用系统时间等随机源初始化随机数生成器
    return random.uniform(0.9990, 1.0001)

def GetRandomNumber():
    random.seed()  # 使用系统时间等随机源初始化随机数生成器
    return random.uniform(-0.1, 0.1)

plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用 SimHei 字体
plt.rcParams['font.weight'] = 'bold'  # 设置字体为粗体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

# 定义季节函数
def season_factor(t):
    # 将一年分为四个季节，每个季节分别为3个月
    season = (t % 12) // 3  # 每个季节3个月
    return season

# Define parameters that change with time
def r(t):
    season = season_factor(t)

    if season == 1:
        increase_rate = 0.2
    elif season == 2:
        increase_rate = 0.24
    elif season == 3:
        increase_rate = 0.25
    else:
        increase_rate = 0.16
    return 2.3 + increase_rate + 2 * np.sin(0.005 * t)

# 定义微分方程组
def evolution_forests_to_farmland(t, y, params):
    P, I, V, B, Bee = y  # 设置初始种群数量
    K, yp_B, a_P, b, I_M, et, yp_I, yp, ks, V_M, mu_V, Ba, B_M, mu_B, r_Bee, Bee_M, sat_Bee, mu_Bee = params
    dPdt = r(t) * P * (1 - P / K) + yp_B * V * P - a_P * P * I  # 植物增长
    dIdt = b * I * (1 - I / I_M) - yp * B * I - et * V * I - yp_I * I + (P / 300) * 15  # 昆虫增长
    dVdt = ks * V * (1 - V / V_M) + et * I * V - mu_V * V  # 蝙蝠增长
    dBdt = Ba * B * (1 - B / B_M) + yp * B * I - mu_B * B  # 鸟类增长
    dBeedt = r_Bee * Bee * (1 - Bee / Bee_M) + sat_Bee * P - mu_Bee * Bee  # 蜜蜂增长
    return [dPdt, dIdt, dVdt, dBdt, dBeedt]

# 初始种群数量
P0 = 450  # 植物
I0 = 120  # 昆虫
V0 = 30  # 蝙蝠
B0 = 45  # 鸟类
Bee0 = 44  # 蜜蜂
y0 = [P0, I0, V0, B0, Bee0]

# 时间点
t = np.linspace(0, 120, 1000)

# 参数设计
params = (
    600,  # K: 植物环境容纳量
    0.007,  # yp_B: 蝙蝠给植物授粉率
    0.005,  # a_P: 植物对昆虫的抗性系数
    1.2,  # b: 昆虫增长率
    200,  # I_M: 昆虫种群最大承载量
    0.003,  # et: 蝙蝠捕虫率
    0.35,  # yp_I: 昆虫自然死亡率
    0.00055,  # yp: 鸟捕虫的效率
    0.25,  # ks: 蝙蝠生长率
    50,  # V_M: 蝙蝠环境容纳量
    0.40,  # mu_V: 蝙蝠自然死亡率
    0.08,  # Ba: 鸟类增长率
    50,  # B_M: 鸟类最大承载量
    0.07,  # mu_B: 鸟自然死亡率
    0.1,  # r_Bee: 蜜蜂生长率
    50,  # Bee_M: 蜜蜂环境容纳量
    0.01,  # sat_Bee: 蜜蜂和植物间的营养关系
    0.08  # mu_Bee: 蜜蜂自然死亡率
)

# 使用 solve_ivp 求解微分方程组，尝试使用 'BDF' 方法
sol = solve_ivp(evolution_forests_to_farmland, [t[0], t[-1]], y0, args=(params,), t_eval=t, method='BDF', atol=1e-6, rtol=1e-3)

# 检查求解结果
if not sol.success:
    print("求解失败:", sol.message)

# 输出维度信息进行调试
print("时间点数量:", len(t))
print("解的维度:", sol.y.shape)

# 绘制结果
plt.figure(figsize=(10, 8))
plt.plot(t, sol.y[0] / P0, label='Current biomass of plants / Initial biomass (P)', color='blue')
plt.plot(t, (sol.y[1] / I0) * GetRandomFactor() + GetRandomNumber() * 0.55 * GetRandomFactor() * np.sin(t), label='Current biomass of insects / Initial biomass (I)', color='orange')
# plt.plot(t, sol.y[2] / V0 + GetRandomNumber() * 0.25 * GetRandomFactor() * np.sin(t), label='Current biomass of bats / Initial biomass (V)', color='red')
plt.plot(t, sol.y[3] / B0 + GetRandomNumber() * 0.1 * GetRandomFactor() * np.sin(t), label='Current biomass of birds / Initial biomass (B)', color='green')
plt.plot(t, sol.y[4] / Bee0+ GetRandomNumber() * 0.15 * GetRandomFactor() * np.sin(t), label='Current biomass of bees / Initial biomass (Bee)', color='purple')

plt.legend()

plt.xlabel('Time (months)', fontweight='bold')
plt.ylabel('Initial biomass / Current biomass', fontweight='bold')
plt.title('Population size over time after farmland ecosystem maturation (no chemicals)', fontweight='bold')
# Set y-axis ticks to 0.1 intervals
plt.yticks(np.arange(0.9, 1.5, 0.1))

# Set x-axis ticks to 12 intervals
plt.xticks(np.arange(0, 121, 12))

plt.show()